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Spatial distribution of the emission intensity in a photonic crystal: Self-interference of 

Bloch eigenwaves 

Dmitry N. ChigrirQ 

Physikalisches Institut, Universitat Bonn, Nussallee 12, D-53115 Bonn, Germany 

The analysis of an angular distribution of the emission intensity of a two-level atom (dipole) in a 
photonic crystal reveals an enhancement of the emission rate in some observation directions. Such 
an enhancement is the result of the bunching of many Bloch eigenwaves with different wave vectors 
in the same direction due to the crystal anisotropy. If a spatial distribution of the emission intensity 
is considered, the interference of these eigenwaves should be taken into account. In this paper, the 
far-field emission pattern of a two-level atom is discussed in the framework of the asymptotic analysis 
the classical macroscopic Green function. Numerical example is given for a two-dimensional square 
lattice of air holes in polymer. The relevance of results for experimental observation is discussed. 

PACS numbers: 42.70.Qs; 42. 25. Fx; 42.50.Pq; 81.05.Zx 



I. INTRODUCTION 

Light emission in a photonic crystal has attracted a 
substantial attention both in theoretical^ E|, B^l-jl, Jlj H , 

l?Jrl^ I^\^d.^ilj22l| studies. To a large extend this 
interest is due to potential perspectives of the emission 
modification and control provided by photonic crystals. 
The inhibition of spontaneous emission is possible within 
a spectral range of a complete photonic bandgap [l], Q , 
where linear propagation of light is prohibited in all spa- 
tial directions. An emission enhancement is a result of 
a long interaction time of an emitter and the radiated 
field, when the emitter is coupled to the slow eigenmode 
0) O HI] or to the strongly localized mode of a defect 
state of a photonic crystal @, [H, [26| . 

It is well known that, the spontaneous decay of an 
excited atom strongly depends on the environment [27| . 
Both the emission rate and the emission directionality 
can be affected. In the simplest case of a two-level atom 
placed in an inhomogeneous medium the emission dy- 
namics can be described by the integro-differential equa- 
tion for the upper state occupation probability amplitude 
8. 281 
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is the projected local density of states (PLDOS) [2a.l3Cj|. 
G (r,r',cj) is the classical macroscopic dyadic Green 



function defined by the inhomogeneous wave equation 



e(r) - V x Vx I G (r, r', u) = - 6 e± (r - r') . 



or 
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Here 5 e , (r — r ) is the e-transverse dyadic delta func- 
tion P. |32|. ro, ojq and d are atom location, transition 
frequency and transition dipole moment, respectively, c 
is a speed of light in vacuum. In the case of a general 
linear, non-magnetic, dielectric medium with arbitrary 
three-dimensional (3D) periodic dielectric function, e(r), 
the Green's function can be expressed in the Bloch mode 
basis as EL [HI 
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Here A„k(r) are Bloch eigenwaves characterized by the 
band index n, the wave vector k„ and the eigenfre- 
quencies w n k- Bloch eigenwaves are solutions of the 
homogeneous wave equation and obey the gauge V • 
[e(r)A ra k(r)1 = 0, normalization and completeness con- 
ditions 0, [32| . The asterisk * and ® denote the complex 
conjugate and the outer tensor product, respectively. A 
positive infinitesimal 7 assures causality. 

In the limit of the weak coupling (Markov approxima- 
tion), a coarse-grained description of the atomic motion 
memory effects can be disregarded [28| and Eq. (Q]) 
yields the familiar exponential decay of the excited state 
with a decay rate given by 
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This approximation gives the correct result for emission 
modification in most of the situations considered in the 
present paper. In the same time, special care should be 
taken for frequencies near the photonic band edges or 
other van Hove singularities, where the memory effects 
become significant and Eq. (TT]) should be analyzed in- 
stead of a direct use of Eqs. (|6l7p . 

In Eqs. (P) and (|6I7|) all parameters of a periodic envi- 
ronment relevant for the atomic evolution are contained 
via the classical Green function and its Bloch mode ex- 
pansion ([5]). The emission rate ^ is proportional to 
the density of available Bloch eigenmodes weighed by the 
coupling strength between the atomic dipole moment and 
the corresponding mode. The emitted intensity |J7J| at the 
point r is determined both by the spectral and angular 
emission rate modification and by the interference of the 
Bloch modes at this point. 

In contrast to the total emission rate modification ((6]) , 
angular-resolved emission experiments, e.g. [HIH^, usu- 
ally probe only a small fraction of the solid angle detect- 
ing the emission modification in a particular direction in 
space. To analyzed such angular-resolved emission ex- 
periments, a fractional power emitted per solid angle in 
space can introduced. It has been shown, that the radi- 
ation pattern of the classical dipole in a photonic crystal 



can demonstrate a strong modification with respect to 
the dipole radiation pattern in vacuum [2lL [33[ ■ For ex- 
ample, in a photonic crystal with incomplete bandgap, 
the angular-resolved emission rate is suppressed in the 
direction of the spatial stopband and strongly enhanced 
in the direction of the group velocity, which is station- 
ary with respect to a small variation of the wave vector 
(photon focusing) [33| . Such an enhancement is the result 
of the bunching of many Bloch eigenwaves with different 
wave vectors in the same spatial direction due to the crys- 
tal anisotropy [34], [35[. For a coherent light source, this 
inevitably leads to the interference of the Bloch eigen- 
waves at the detector plane and additional modification 
of the spatial distribution of the emission intensity. 

In this paper, the physical picture of the interference 
fringes formation in the far-field emission pattern of the 
two-level atom placed in a periodic medium is considered. 
The paper is organized as follows. The evaluation of the 
asymptotic form of the emitted intensity |(7j) is given in 
sectionlTTlin the radiation zone. The physical explanation 
and the relevance of results for experimental observation 
are discussed in section Mil In section IIVI a numerical 
example is given for a two-dimensional polymer photonic 
crystal. Section |V] summarized the main results of the 
paper. 

II. ASYMPTOTIC FORM OF EMITTED 
INTENSITY 

By taking into account the Bloch theorem, A n k(r) = 
a„k(r)e lk "' r , where a n k(r) is a lattice periodic function, 
and changing the fc-space summation to the correspond- 
ing fc-space integral, ^ k — > (V/8tt 3 ^ J d 3 k, Eq. ([7]) can 
be expressed as 
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For large |x| = |r — ro| the exponential function in the integral (|5|) oscillates rapidly. To evaluate the integral, the 
method of stationary phase can be used. As it was shown in [33j , the principal contribution to the integral comes from 
the regions of the iso-frequency surface in the wave vector space, at which the eigenwave group velocity is parallel to 
observation direction x. Then, the integral in Eq. ([5]) can be transformed to the form 
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where V^ ik is the group velocity of the Bloch eigenwave, the integration is over the iso-frequency surface w nk — ujq 
and summation is taken over all stationary eigenwaves v with a group velocity vector pointing in the observation 
direction x = r ro- 

The result of the integration in ([9]) depends on the local topology of iso-frequency surface w„k = uiq. It is convenient 
to introduce the local curvilinear coordinates £j with the origin at the iso-frequency surface and with one of the 
coordinate aligned perpendicular to it, e.g., £3. Then, a function h (£1, £2) = k„ • x can be expanded in a series near 
the wave vector of the eigenwave (p, n, k): 
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and x is a unit vector in the observation direction x = r — i"o- If the iso- frequency surface has a non- vanishing Gaussian 
curvature in the vicinity of the wave vector k^, only quadratic terms in the expansion (|10p can be kept, leading to 
the following asymptotic from of the far-field intensity© [33] 
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where — oc^a^ determines the Gaussian curva- 

ture of the iso-frequency surface at the point k„ = kj y t 
(stationary point) and summation is over all stationary 
points with x • V,^ k > 0. 

The emission intensity far from an atom is proportional 
to the inverse Gaussian curvature of the iso-frequency 
surface, ~ > an( l to the inverse square of the 

distance between the source and the observation point, 
~ |x| . The asymptotic energy flux shows the necessary 
amount of decrease with distance (~ |x| 2 ), providing a 
finite value of the energy flux in any finite interval of 
a solid angle, that is, assuming non-vanishing Gaussian 
curvature. A vanishing curvature formally implies an in- 
finite flux along the corresponding observation direction, 
leading to the photon focusing phenomenon [H, [36|, H3] • 

Strictly speaking, the asymptotic behavior of the emis- 
sion intensity is valid only if quadratic terms in the 
expansion (TIT)|) do not vanish, so that all higher order 
terms in the expansion can be neglected. A parabolic 
point of the iso-frequency surface is an example of van- 
ishing quadratic terms in (| 10|) . Generally, the Gaussian 
curvature is zero at a parabolic point and one (both) of 



the principal curvatures of the iso-frequency surface is 
zero. Actually, at parabolic points the asymptotic be- 
havior of the emitted intensity i.e., the dependence of 
the intensity on the inverse distance (~ |x| 1 ) changes 
to the power of the inverse distance. 

As an illustration one can consider the simple parabolic 
point ko = k° in the vicinity of which the function 
M^iifo) nas t ne expansion [38| : 

6) =k -x + ia£ ! + i/3£ 2 3 , a = a° n , = P° ln , 
2 b 

(12) 

where Xo is the unit vector in the direction normal to the 
iso-frequency surface at the parabolic point ko . The local 
curvilinear coordinates £j has the origin at the parabolic 
point ko, with the coordinates £i and £2 aligned along the 
directions of the principal curvatures of the iso-frequency 
surface at this point and with the coordinate £3 aligned 
along Xo. For the parabolic point ko ([12"]) one of the prin- 
cipal curvature vanishes (a^ = 0), while another princi- 
pal curvature remain non-zero. Using the expansion (|12|) 
the asymptotic form of the intensity © is given by: 
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where Ao(r) = ao(r)e lk °' r and Vo are the Bloch mode and the group velocity associated with the parabolic point ko, 
respectively. Calculating integrals in (fT3"|) , leads to 
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for direction £2- Here T (|) is the Gamma function. Now, combining (|14p and (| 1 5|) the following expression for the 
asymptotic vector potential associated with the parabolic point (|12p can be obtained [H, H(| : 
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The emission intensity associated with a parabolic to the usual inverse square law |x| for other directions. 



point falls off with the distance as |x| 



-5/3 



in contrast 
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If there are no additional singularities on the parabolic 

1 li 

line, the product \a\ \(3\ ' ~ K, where K is the Gaus- 
sian curvature at an arbitrary point of the iso-frequency 
surface. Then the emission intensity is proportional to 
K^ 1 \f3\~ x ^ . Thus, at large |x|, the energy flux along 
the direction corresponding to a parabolic point on the 
iso-frequency surface exceeds the energy flux along the 
direction corresponding to an elliptical point in the ratio 

|x| 1/3 |^| 1/6 . 

The expression {jTHJ) gives the asymptotic emitted in- 
tensity in the direction xo associated with a parabolic 
point on the iso-frequency surface ko , so in the direction 
of the group velocity at the parabolic point. Now, the 
asymptotic intensity for directions x near the direction 
of that group velocity will be calculated. As before, the 
origin of the coordinates & is chosen at the parabolic 



point, where the direction of observation x is coincides 
with direction x . It is assumed, that the principal cur- 
vature vanishes in the £2 direction. Let the position x 
(x||x) be described by coordinates Xj. Then, since x is 
nearly parallel to Xo one have from (| 1 2|) [38] : 

6 = \< + 

and 

k n -x«k -x + £ixi + fa 2 +(\a£l + \p§\ X3 . (17) 

Using expansion <| 1 7|) the asymptotic form of the intensity 
© is given by 
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where the fact that |x| and X3 are approximately equal was used. The integral in (|18p over £1 is calculated simply to 
be 



/ ? 1 1 1/2 1 ii 

J-00 \a\ 1 |x| 

while the integral over £2 results in 
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where Ai is the Airy function. Then the asymptotic emitted intensity (|18p is finally given by 
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As in the case of the asymptotic intensity (fT6|) . en- 
ergy flux in the direction x falls off as |x|~ 5 ^ 3 and ex- 
ceeds the energy flux along the other directions in the 
ratio |x| 1 ^ 3 l/JI 1 ^ 6 . The dependence of the energy flux in 
the plane of the vanishing principal curvature is given 
by the square of the Airy function [Ai (a^/a)] 2 - When 
x-ija is positive, the energy flux is small and exponen- 
tially drops while the angle between direction of obser- 
vation and direction corresponded to the parabolic point 
increases (Fig. [IJ. For negative x^ja, the flux oscillates 
rapidly and has a mean value averaged over one cycle pro- 
portional to ~ (x2/a)^ 1 ^ 2 (FigQ]). A mean value of the 
energy flux is then proportional to |a| _1 |/3| -1 ^ 2 ~ K^ 1 
and |x|~ 2 and coincides with the asymptotic energy flux 



associated with an elliptical point of iso-frequency sur- 
face (TTT|) . demonstrating focusing of the energy flux in 
the direction corresponding to the parabolic point of the 
iso-frequency surface. 



III. INTERFERENCE OF BLOCH 
EIGENWAVES 

The vanishing curvature of the iso-frequency surface 
results in the folds of the wave front (wave surface) [H[ . 
Then, for the direction near the fold of the wave surface 
the field is a superposition of several Bloch eigenwaves 
(Fig. [5J). In the far-field, where the source-to-detector 
distance is much larger than the source size and the wave- 
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Figure 1: Plot of Airy function [Ai (0:2/0)] 2 (solid line) and 
its mean value (l/47r) _1 (x2/a)~ 1 ^ 2 (dashed line). The dashed 
line can be considered as the "geometrical optics" approxima- 
tion of the Airy function. 



length, the part of the wave front limited to the small 
solid angle can be approximated as a Bloch eigenwave 
with the group velocity within this angle. If there is a 
relative difference in the lengths or directions of the wave 
vectors of Bloch eigenwaves, the eigenwaves can interfere 
yielding an oscillations in the energy flux distribution. 
This can be already seen from the general expression for 
the emitted intensity (0. 

Two general conditions are required for the interfer- 
ence to occur. The polarization states of the Bloch 
eigenwaves must be nonorthogonal and the Bloch eigen- 
waves must overlap in space (4l| . This kind of inter- 
ference of the Bloch eigenwaves will be called further a 
self-interference, to stress that the field produced by the 
light source inside a photonic crystal can interfere with 
itself producing an interference patten in the energy flux 
distribution. A similar self-interference effect also hap- 
pens in the case of ballistic phononspropagation in an 
acoustically anisotropic crystals [13, E2] • 

For a more qualitative measure of the self-interference 
effect the iso-frequency surface superimposed on a pho- 
tonic crystal is considered in figure [3J Here evaluation 
presented by Hauser et al. [42j for the self-interference of 
ultrasound in a crystal is followed. In the figure [3l dots 
are parabolic points of zero curvature. The light source 
is located near the bottom surface of the photonic crystal 
and generates uniform distribution of wave vectors. The 
Bloch eigenwave with wave vector ko propagates with the 
group velocity V , normal to the iso-frequency surface at 
ko, arriving at the point Ro on the opposite surface of the 
crystal. Near the parabolic point the iso-frequency sur- 
face are practically flat, neighboring wave vectors have 
nearly the same group velocity. This gives rise to the 
high-intensity caustic in the detected intensity distribu- 
tion (Fig. [3]- left). If the detector is moved to a point 
Ri slightly away from Ro, two distinct Bloch eigenwaves 
with different wave vectors k' x and k' 2 near ko arrive at 
the detector (Fig. [3J- right). If the surface were perfectly 
flat near ko, then k^-Ri = k 2 -Ri, and the two eigenwaves 
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Angle (degree) 

Figure 2: Diagram illustrating the self-interference of the 
Bloch eigenwaves in a photonic crystal. A section of the wave 
surface with a fold is presented (thick solid line). The asymp- 
totic intensity (solid line) displays oscillations as the angle 
passes the fold section of the wave surface. The asymptotic 
intensity in the "geometrical optics" limit demonstrates focus- 
ing caustics in the direction of the folds (dashed line) . Within 
a small solid angle the far-field of a point source consists on 
superposition of three Bloch eigenwaves with different wave 
vectors. These Bloch eigenwaves interfere leading to oscilla- 
tions in the intensity distribution. Only two wave vectors are 
illustrated for clarity. 



would always remain in phase at the detector, interfer- 
ing constructively. In reality the iso-frequency surface 
is curved near the parabolic point, so as Ri is rotated 
downward the corresponding waves begin to interfere de- 
structively, producing an Airy pattern (Figs. [T][3]). If k^ 
and k 2 are close two ko , and that q = k^ — ko ~ k 2 — ko , 
then destructive interference will take place, if the to- 
tal phase difference of the light as it travels through the 
sample, 2q • R 1; is an odd integer multiple of tt. 

Strictly speaking, there is one more Bloch eigenwave 
following in the observation direction in the fold region 
of the wave surface (Figs. OE])- This eigenwave is de- 
picted as Vi and V3 in the left and right panels of 
figure [3l respectively. To obtain a complete picture of 
the self-interference near the fold of the wave surface, 
a three-wave interference should be taken into account, 
which would lead to more complicated interference pat- 
terns in the intensity distribution. Here, the influence 
of this third eigenwave is neglected for simplicity, hence 
this eigenwave usually has a relatively small group veloc- 
ity and could arrive at the detector too late to interfere 
with the eigenwaves k^ and k 2 . 

From the perspective of the self-interference effect, the 
mean value of the asymptotic energy flux (|19[) averaged 
over one cycle of Airy oscillations can be viewed as a 
"geometrical optics" approximation of the actual energy 
flux. This approximation then corresponds to the ray 
description of wave propagation, where the energy flux 
is simply proportional to the density of rays crossing a 
detector surface. In this picture the interference among 
different rays is neglected. Then, the emission rate en- 
hancement in the focusing direction can also be inter- 
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Figure 3: Schematic of self-interference in a photonic crys- 
tal. The iso-frequency surface is a sketch of the iso-frequency 
surface of a real 2D photonic crystal. The plot on the top is 
the detected intensity distribution. Left: For the wave vec- 
tor ko, the group velocity is Vo and the wave arrives at Ro. 
Right: There are two eigenwaves with different wave vectors 
kj, which group velocity points in the same direction Ri. The 
difference in the wave vectors of these eigenwaves results in 
the waves arriving at Ri with different phases and leads to 
interference. 

preted as the relative increase of the rays density or as an 
increased probability of the photon emission in this ob- 
servation direction. As it has been mentioned above, the 
mean value of the asymptotic flux (fl9|) coincides with the 
asymptotic energy flux (|11|) derived for elliptical points 
of the iso-frequency surface. So, the asymptotic energy 
flux (jTTJ) corresponded to an elliptical point can be also 
considered as a "geometrical optics" approximation, and 
can be used for all points of the iso-frequency surface 
within this approximation. 

In a typical experiment the differences between the en- 
ergy flux (|19p and its "geometrical optics" approximation 
(fTTj) will be reduced by the effect of the finite size of the 
light source and the detector. It is clear that if the lin- 
ear dimensions of the source area and detector are L, 
the intensity is averaged over xi values with a spread of 
L. To see the oscillations of the energy flux one therefore 
needs L < A8 x Ri , where R\ is the distance between the 
source and the detector and A8 is an angular separation 
of the fringes of intensity distribution. 

To estimate this angular separation, Bloch eigenwaves 
k'j and are further approximated by plane waves. 
Then their superposition at the detector position Ri , as- 
suming that they have the same polarization, is (42J: 

e iki Rl + e <kJ = 2 cos (Ak • R x ) e ik ° Rl , 

which is a plane wave with average wave vector ko = 
(k', + k^) /2, modulated by a cosine function with effec- 
tive wave vector Ak/2 = (k' 2 - ki) /2. When Ak Ri = 
AkttRi = 7r, the waves interfere destructively at the de- 
tector. To estimate Afen a local Cartesian coordinate 
system £j with the origin at k and £ 3 along V is chosen 
as it is shown in the figure[3]-left. Then, the iso-frequency 
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Figure 4: Photonic band structure of the square lattice of air 
holes made in a polymer. Polymer has the refractive index 
f.56. Radius of holes is r = 0.15d, where d is the lattice 
period. The band structure is given for TM polarization. The 
frequency is normalized to fl = uid/2-nc = d/X. c is the speed 
of light in the vacuum. Insets show the first Brillouin zone 
(right) and a part of the lattice (left). 



surface near the parabolic point can by parametrized as 
6 = -a£!/ fc o and Afe|| = -2£ 3 = 2aff//c 2 0. There- 
fore, the first minimum in the intensity will occur when 

6 = K/2 a i?!) 1/3 ^i?r 1/3 A- 2 / 3 , 

where A = 2ir/ko is the average wavelength. Finally, the 
coordinate-space angle between the intensity maximum 
and the first minimum is given by [42| : 

A9 » |VJ - V | /V = 3a (7r/2ak R 1 ) 2/3 ~ (A/i?r) 2/3 . 

Then for optical wavelength, e.g., 500 nm, and a dis- 
tance to the detector of 1 cm the linear dimension of 
the light source and the spatial resolution of the detector 
should be smaller than 10 fim. So, in most experiments 
the "geometrical optics" approximation pip reasonably 
represents an asymptotic emission intensity of the light 
source inside a photonic crystal. 



IV. NUMERICAL EXAMPLE 

In this section different approximations of the light 
emission pattern, discussed in section [ill are compared. 
Numerical calculations are done for a point source placed 
inside a two-dimensional polymer photonic crystal. A 
point source produces an isotropic and uniform distribu- 
tion of wave vectors k n with the frequency ujq. Then, the 
asymptotic field in ([9]) and (jTTJ) should be averaged over 
the dipole moment orientation, which yields a factor of 
|d|/3. 

As it was pointed out in section [TlJ the main contribu- 
tion to the far-field of a point source inside a photonic 
crystal comes from the vicinity of the wave vector of the 
eigenmodes with the group velocity in the observation 
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tonic crystal (Fig. [4| for the normalized frequency Q = 0.333. 

Parabolic point are marked by the black dots. The first Bril- Figure 6: Wave contour corresponded to the normalized fre- 
louin zone of the lattice is plotted in order to show the spatial quency Q = 0.333. The group velocity is plotted in the units 
relation between zone boundary and iso-frequency contours. of the speed of light in vacuum. The directions corresponded 

to the folds of the wave contour are shown. 



direction. That means that the far-field emission inten- 
sity of a point source is mainly given by the square of the 
integral in © 



2l, e ik„(r-r ; 



d 2 k. 



(21) 



In what follows, the contribution only from one photonic 
band is considered. In the "geometrical optics" approxi- 
mation the main contribution to the far-field emis- 
sion intensity is then given by an inverse Gaussian cur- 
vature of the iso-frequency surface, 



(22) 



The angular distribution of the far-field intensity de- 
pends on the topology of the iso-frequency surface of 
the crystal at the emission frequency (|21ti22|) . In what 
follows, an infinite two-dimensional square lattice of air 
holes in a polymer background is considered. Polymer 
has the refractive index n — 1.56, radius of holes is 
r = 0.l5d, where d is the lattice period. The consid- 
eration is limited to the in-plane propagation of the TM 
mode of the crystal. The photonic band structure of such 
a photonic crystal is presented in the figure [4j The band 
structure has been calculated using the plane wave ex- 
pansion method (43j . 

The iso-frequency contours for the normalized fre- 
quency £1 = 0.333 is presented in figure[5] The frequency 
belongs to the first photonic band and it is within the 
first stopband in the TX direction of the crystal. To plot 
an iso-frequency contour, the photonic band structure for 
all wave vectors within the irreducible Brillouin zone was 
calculated and then an equation w(k) = loq was solved 
for a given frequency ujq. The iso-frequency contour is an 
open contour and has alternating regions of the Gaussian 




Figure 7: Angular distribution of radiative power correspond- 
ing to the normalized frequency Q = 0.333. The directions of 
infinite radiative power (caustic) coincide with the directions 
of the folds of the wave contour (Fig. [6j . 



curvature with a different sign. Parabolic points, where 
the Gaussian curvature vanishes, are marked by black 
dots in the figure [5] The vanishing curvature results in 
the folds of the wave contour and in the focusing of the 
light in the folds direction (33|. The wave contour corre- 
sponding to the iso-frequency Q — 0.333 is presented in 
the figure O A pair of the parabolic points in the first 
quarter of the Brillouin zone results in a cuspidal struc- 
ture of the wave contours in the first quarter of the coor- 
dinate space. In the figure [7] the polar plot of the main 
contribution to the far-field intensity in the "geometrical 
optics" approximation (|22|) is presented. The energy flux 
is strongly anisotropic, showing relatively small intensity 
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Figure 8: "Geometrical optics" (|22p and "wave optics" (12ip 
approximations of the far-field emission intensity. Normalized 
frequency is Q = 0.333. The distance between a point source 
and a detector is 100 lattice periods. The top panel is for the 
"geometrical optics" approximation. The middle and bottom 
panels are for the "wave optics" approximations. See text for 
more details. 



in the directions of the stopband, and infinite intensity 
(caustics) in the directions of the folds. 

In figure [S] the comparison of the "geometrical op- 
tics" (f2"2"j) and "wave optics" (f2"Tj) approximations of the 
far-field emission intensity is given. A normalized in- 
verse Gaussian curvature of the iso-frequency surface 
O = 0.333 is presented in the top panel. Focusing di- 
rections, 22° apart from the [10] direction of the square 
lattice, are clearly seen. The integral (|2"Tj) evaluated for 
the distance 100 period apart from the point source is 
given in the middle and the bottom panels of figure [5] 
The normalized intensity distribution presented in the 
middle panel was calculated by reducing the integration 
limits in (|2"Tj) to the close neighborhood of the parabolic 
point of the iso-frequency surface. Then the result of the 
integration is similar to one in (|19p and an angular dis- 




0.0 0.25 0.5 

Figure 9: An asymptotic map of the intensity distribution 
inside a 50 x 50 photonic crystal. Normalized frequency is 
fl — 0.333. The structure of the crystal is superimposed on 
the field map. Folds directions are shown by black lines. 




0.0 0.005 0.01 

Figure 10: FDTD calculation. Map of the modulus of the 
Poynting vector field for a 50 x 50 rod photonic crystal excited 
by a point isotropic source with the normalized frequency 
Q — 0.333. The asymptotic directions of photon focusing 
caustics are shown as black lines. 



tribution of the emission intensity resembles the square 
of the Airy function. Actually, this approximation takes 
into account an interference of only two Bloch eigenwaves 
in the fold region. If the three wave interference is taken 
into account, by extending the integration limits in (I21|) 
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Figure 11: The comparison of the intensity distribution for 
the asymptotic (solid line) and the FDTD (dashed line) cal- 
culations. The distance between the point source and the 
detector is 20 lattice periods. 



over all iso-frequency surface in the first quarter of the 
Brillouin zone, a more complex interference pattern ap- 
pears in the angular emission intensity distribution (Fig. 
[5]-bottom). Both "wave optics" approximations show an 
intensity enhancement along fold directions. 

In figure [H] a two dimensional map of the intensity 
distribution inside a 50 x 50 photonic crystal is pre- 
sented. The intensity distribution was calculated using 
Eq. (j2~Tj) by integration over complete iso-frequency con- 
tour £7 = 0.333. The structure of the crystal is superim- 
posed on the field map. Folds directions are shown by 
black lines. The focusing of the light in the fold direc- 
tion together with an Airy-like oscillations between folds 
directions are clearly seen in figure [5] 

To substantiate an asymptotic analysis, the finite 
difference time domain (FDTD) calculations was done 
[USE!. The simulated structure was a 50 x 50 lattice. 
The crystal is surrounded by an extra 2d wide layer of 
polymer. The simulation domain was discretized into 
squares with a side A = d/32. The total simulation 
region was 1728 x 1728 cells plus 8-cell wide perfectly 
matched layer (PML) [46]. The point isotropic light 



source was modeled by a current density source [44 145 1 
with a homogeneous spacial dependence and sinusoidal 
temporal dependence of the signal. 

In figure [TU] the map of the modulus of the Poynt- 
ing vector field is shown, when the crystal is excited by 
a point isotropic source. The point source is placed in 
the middle of the crystal. A field map is shown for one 



instant time step. The snap-shots were captured after 
10000 time steps, where the time step was 4.38 x 10~ 17 
s (0.99 of the Courant value). The structure of the crys- 
tal is superimposed on the field map. One can see, that 
the emitted light is focused in the directions of the folds 
(black lines). Moreover, an interference pattern between 
the folds directions is in a reasonable agreement with 
the interference pattern predicted using the asymptotic 
analysis (Fig. [9]) . For the FDTD calculations, a periodic 
modulation of the intensity in the radial direction will 
go away if time averaging is performed. The comparison 
of the intensity distribution, 20 periods apart from the 
point source is given in figure [TT] A reasonable agreement 
between interference minima and maxima positions for 
the asymptotic (solid line) and the FDTD (dashed line) 
calculations is shown. The disagreement in the absolute 
values of the angular intensity distributions is mainly be- 
cause of the prefactor of the integral in which was 
neglected in (|2"Tj) for simplicity. 



V. SUMMARY 

It was shown, that the intensity modulation of the an- 
gular resolved emission spectra is not only due to the 
emission rate modification, but also is the result of the 
interference of several photonic crystal eigenmodes with 
different wave vectors approaching detector at the same 
moment of time. Using an asymptotic analysis of classi- 
cal Green function, "geometrical optics" and "wave op- 
tics" approximations of the emitted intensity due to a 
two-level atom were introduced in the radiation zone. 
The physical reasons for the interference pattern forma- 
tion and the possibilities of experimental observation of 
them were discussed. A numerical example was given 
in the case of polymer two-dimensional photonic crystal. 
It was shown that rigorous FDTD calculations arc in a 
reasonable agreement with the developed approximate 
analysis. 



Acknowledgments 

This work was partially supported by the EU-IST 
project APPTech IST-2000-29321, the BMBF project 
PCOC 01 BK 253 and the DFG Research Unit 557. 



[1] V. P. Bykov, Soviet Physics - JETP 35, 269 (1972). 
[2] E. Yablonovitch, Phys. Rev. Lett. 58, 2059 (1987). 
[3] S. John and J. Wang, Phys. Rev. Lett. 64, 2418 (1990). 
[4] J. P. Dowling and C. M. Bowden, Phys. Rev. A 46, 612 
(1992). 

[5] T. Suzuki and P. K. L. Yu, J. Opt. Soc. Am. B 12, 570 
(1995). 

[6] K. Sakoda and K. Ohtaka, Phys. Rev. B 54, 5732 (1996). 
[7] S. Nojima, Jpn. J. Appl. Phys. 2 37, L565 (1998). 



[8] K. Busch, N. Vats, S. John, and B. C. Sanders, Phys. 

Rev. E 62, 4251 (2000). 
[9] Y. Xu, R. K. Lee, and A. Yariv, Phys. Rev. A 61, 033807 

(2000). 

[10] Z. Y. Li, L. L. Lin, and Z. Q. Zhang, Phys. Rev. Lett. 

84, 4341 (2000). 
[11] V. Lousse, J. P. Vigneron, X. Bouju, and J. M. 

Vigoureux, Phys. Rev. B 64, 201104R (2001). 
[12] C. Hermann and O. Hess, J. Opt. Soc. Am. B 19, 3013 



10 



(2002). 

[13] J. P. Dowling, M. Scalora, M. J. Bloemer, and C. M. 

Bowden, J. Appl. Phys. 75, 1896 (1994). 
[14] H. Hirayama, T. Hamano, and Y. Aoyagi, Appl. Phys. 

Lett 69, 791 (1996). 
[15] M. Boroditsky, R. Vrijen, T. F Krauss, R. Coccioli, 

R. Bhat, and E. Yablonovitch, IEEE J. Lightwave Tech- 

nol. 17, 2096 (1999). 
[16] S. V. Gaponenko, V. N. Bogomolov, E. P. Petrov, A. 

M. Kapitonov, A. A. Eychmueller, A. L. Rogach, I. I. 

Kalosha, F. Gindele, and U. Woggon, J. Lightwave Tech- 

nol. 17, 2128 (1999). 
[17] B. Temelkuran, M. Bayindir, E. Ozbay, R. Biswas, M. 

M. Sigalas, G. Tuttle, and K. M. Ho, J. Appl. Phys 87, 

603 (2002). 

[18] H. P. Schriemer, H. M. van Driel, A. F. Koenderink, and 
W. L. Vos, Phys. Rev. A 63, 011801R (2000). 

[19] S. G. Romanov, T. Maka, C. M. Sotomayor Torres, 
M. Muller, and R. Zentel, Appl. Phys. Lett. 79, 731 

(2001) . 

[20] A. F. Koenderink, L Bechger, H. P. Schriemer, A. La- 
gendijk, and W. L. Vos, Phys. Rev. Lett. 88, 143903 

(2002) . 

[21] S. G. Romanov, D. N. Chigrin, V. G. Solovyev, T. Maka, 
N. Gaponik, A. Eychmuller, A. L. Rogach, and So- 
tomayor Torres, Phys. Stat. Sol. 197, 662 (2003). 

[22] A. F. Koenderink, L Bechger, A. Lagendijk, and W. L. 
Vos, Phys. Stat. Sol. 197, 648 (2003). 

[23] K. Sakoda, Opt. Express 4, 167 (1999). 

[24] J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Pho- 
tonic crystals: molding the flow of light (Princeton Uni- 
versity Press, Princeton NJ, 1995). 

[25] R. Coccioloi, M. Boroditsky, K. W. Kim, Y. Rahmat- 
Samii, and E. Yablonovitch, Inst. Elct. Eng. Pros.- 
Optoclectron. 145, 391 (1998). 

[26] O. J. Painter, A. Husain, A. Schercr, J. D. O'Brien, 
I. Kim, and P. D. Dapkus, J. Lightwave Tech. 17, 2081 
(1999). 



[27] E. M. Purcell, Phys. Rev. 69, 681 (1946). 

[28] H. T. Dung, L. Knoll, and D.-G. Welsch, Phys. Rev. A 

62, 053804 (2000). 
[29] R. Sprik, B. A. Van Tiggelen, and A. Lagendijk, Euro- 

phys. Lett. 35, 265 (1996). 
[30] R. C. McPhedran, L. C. Botten, J. McOrist, A. A. 

Asatryan, C. M. de Sterke, and N. A. Nicorovici, Phys. 

Rev. E 69, 016609 (2004). 
[31] K. Sakoda, Optical Properties of Photonic Crystals 

(Springer, Berlin, 2001). 
[32] R. J. Glauber and M. Lewenstein, Phys. Rev. A 43, 467 

(1991). 

[33] D. N. Chigrin, Phys. Rev. E 70, 056611 (2004). 

[34] R. Zengerle, J. Mod. Optics 34, 1589 (1987). 

[35] P. St. J. Russell, Appl. Phys. B: Photophysics & Laser 

Chemistry B39, 231 (1986). 
[36] P. Etchegoin and R. T. Phillips, Phys. Rev. B 53, 12674 

(1996). 

[37] D. N. Chigrin and C. M. Sotomayor Torres, Optics and 

Spectroscopy 91, 484 (2001). 
[38] A. M. Kossevich, The Crystal Lattice (Wiley, Berlin, 

1999). 

[39] A. A. Maradudin, in Phonons and Phonon Interactions, 
edited by T. A. Bak (Benjamin, New York, 1964), pp. 
424-504. 

[40] H. J. Maris, Phys. Rev. B 28, 7033 (1983). 

[41] P. St. J. Russell, Phys. Rev. A 33, 3232 (1986). 

[42] M. R. Hauser, R. L. Weaver, and J. P. Wolfe, Phys. Rev. 

A 33, 3232 (1986). 
[43] S. G. Johnson and J. D. Joannopoulos, Opt. Express 8, 

173 (2001). 

[44] A. Taflove, Computational Electrodynamics: The Finite- 
Difference Time-Domain Method (Artech House, Nor- 
wood, 1995). 

[45] D. M. Sullivan, Electromagnetic Simulation Using the 

FDTD Method (IEEE Press, New York, 2002). 
[46] J. P. Berenger, J. Comp. Phys. 114, 185 (1994). 



